A Supersymmetric Approach to Excited States via Quantum Monte Carlo 
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We present here a supersymmetric (SUSY) approach for determining excitation energies within the 
context of a quantum Monte Carlo scheme. By using the fact that SUSY quantum mechanics gives 
rises to a series of isospectral Hamiltonians, we show that Monte Carlo ground-state calculations 
in the SUSY partners can be used to reconstruct accurately both the spectrum and states of an 
arbitrary Schrodinger equation. Since the ground-state of each partner potential is node-less, we 
avoid any "node" -problem typically associated with the Monte Carlo technique. While we provide 
an example of using this approach to determine the tunneling states in a double-well potential, the 
method is applicable to any ID potential problem. We conclude by discussing the extension to 
higher dimensions. 



I. INTRODUCTION 

The variational Monte Carlo (VMC) technique is a 
powerful way to estimate the ground state of a quantum 
mechanical system. The basic idea is that one can use 
the variational principle to minimize energy expectation 
value with respect to a set of parameters {a} 



E(a) 



J \^p{x,a)\'^{Hij)/'il;{x,a))dx 
J \ip{x, a)\^dx 



(1) 



Following the Monte Carlo method for evaluating inte- 
grals, one intreprets 



p{x)dx 



\ip{x, a)\'^dx 
J \i^{x, a)\'^dx 



(2) 



as a probability distribution function. Typically, one 
assumes a functional form for the trial wave function, 
ip{x,a) and the numerical advantage is that one can eval- 
uate the energy integral by simply evaluating ijj{x,a). 
The method becomes variational when one then adjusts 
the parameters to optimize the trial wave function. Since 
the spectrum of H is bounded from below, the opti- 
mized trial wave-function provides a best approximation 
to the true ground state of the system. However, since 
p{x) = lijjix, a)^ is a positive definite function, this pro- 
cedure fails if the system has nodes or if the position of 
the nodes is determined by the parameters. One can in 
principle obtain excitation energies by constraining the 
trial function to have a fixed set of nodes perhaps deter- 
mined by symmetry. 

Given that (VMC) is a robust technique for ground 
states, it would be highly desirable if the technique could 
be extended to facilitate the calculation of excited states. 
In this paper, we present such an extension (albeit in 
one dimension) using supersymmetric (SUSY) quantum 
mechanics. The underlying mathematical idea behind 



SUSY is that every Hamiltonian Hi — T + Vi has a 
partner Hamiltonian, H2 = T + V2 (T being the kinetic 
energy operator) in which the spectrum of Hi and H2 are 
identical for all states above the ground state of Hi . That 
is to say, the ground state of H2 has the same energy as 
the first excited state of Hi and so on. This hierarchy 
of related Hamiltonians and the algebra associated with 
the SUSY operators present a powerful formal approach 
to determine the energy spectra for a wide number of 
systems. [H, 0, i, i, ill, 0, i i, M M To date, little 
has been done exploiting SUSY as a way to develop new 
numerical techniques. 

In this paper, we shall use the ideas of SUSY-QM to 
develop a Monte Carlo-like scheme for computing the 
tunneling splittings in a symmetric double-well poten- 
tial. While the model can be solved solved using other 
techniques, this provides a useful proof of principle for 
our approach. We find that the the SUSY/VMC combi- 
nation provides a useful and accurate way to obtain the 
tunneling splitting for this system and excited state wave 
function for this system. While our current focus is on 
a one-dimensional system, we conclude by commenting 
upon how the technique can be extended to multi-particle 
systems and to higher dimension. In short, our results 
strongly suggest that this approach can be brought to 
bear on a more general class of problems involving mul- 
tiple degrees of freedom. Surprisingly, the connection 
between the Monte Carlo technique and the SUSY hier- 
archy has not been exploited until this paper. 



II. SUPERSYMMETRIC QUANTUM 
MECHANICS 

Supersymmetric quantum mechanics (SUSY-QM) is 
obtained by factoring the Schrodinger equation into the 
form [il,[l|[i3 
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Hi; = AUv4^) = 
using the operators 

h 



A 



--d^ + W 



(3) 
(4) 



2 



/2m 



=5, + W. 



(5) 



bmce we can impose At}}),'' = 0, we can immediately 
write that 



W{x) = -- 



'2rn 



dx In i/'o- 



(6) 



W{x) is the superpotential which is related to the physical 
potential by a Riccati equation. 



V{x) = M^2(a;) 



'2m 



W'{x). 



(7) 



The SUSY factorization of the Schrodinger can always 
be applied in one-dimension. 

From this point on we label the original Hamiltonian 
operator and its associated potential, states, and energies 
as Hi, Vi, ipn'' and ErP . One can also define a partner 
Hamiltonian, H2 = AA'' with a corresponding potential 



V2 = 



/2m 



W'{x). 



(8) 



All of this seems rather circular until one recognizes that 
Vi and its partner potential, V2, give rise to a common 
set of energy eigenvalues. This principle result of SUSY 
can be seen by first considering an arbitrary stationary 
solution of Hi, 



This implies that (Aipl^^) is an eigenstate of H2 with 
energy En smce 

H^iA^P^P) = AAtA^(i) = (10) 

Likewise, the Schrodinger equation involving the partner 



potential H21P 



(2) 



(2) (2) 

En ipn implies that 



AtAAt^(f) . i/i(At^(f)) = (At^W). (11) 



This (along with Eo^' — ) allows one to conclude that 
the eigenenergies and eigenf unctions of Hi and H2 are 
related in the following way: E^^^ = En\ 



(2) _ 



E. 



(1) 

n+l 



--Aipl^li, and V'i+i 



E^, 



(2) 



=Atvi2) (12) 



for n > 0. [32] Thus, the ground state of H2 has the same 

(2) 

energy as the first excited state of Hi . If this state ipo 
is assumed to be node-less, then ip[^'' cx A^'ip'^'' will have 
a single node. We can repeat this analysis and show that 
H2 is partnered with another Hamiltonian, H3 whose 
ground state is isoenergetic with the first excited state 
of H2 and thus isoenergetic with the second excited state 
of the original Hi. This hierarchy of partners persists 
until all of the bound states of Hi are exhausted. 



III. ADAPTIVE MONTE CARLO 

Having defined the basic terms of SUSY quantum me- 
chanics, let us presume that one can determine an accu- 
rate approximation to the ground state density po^^ (x) of 
Hamiltonian Hi . One can then use this to determine the 
superpotential using the Riccati transform 



(1) _ 



/2m dx 



and the partner potential 



V2^Vi 



d In 
2m dx^ 



(1) 



(13) 



(14) 



Certainly, our ability to compute the energy of the 
ground state of the partner potential V2 depends on hav- 
ing first obtained an accurate estimate of the ground- 
state density associated with the original Vi . 

For this we turn to an adaptive Variational Monte 
Carlo approach developed by Maddox and Bittner.fTsj 
Here, we assume we can write the trial density as a sum 
over N Gaussian approximate functions 



Pt{x) = y^G'„(x, c„). 



(15) 



parameterized by their amplitude, center, and width. 



-C„2(2 



n3) 



(16) 



(9) This trial density then is used to compute the energy 



E[pt] = Tr[pTE{x)] = {Vi) + {Q[pt]) 
where (5[pt] is the Bohm quantum potential, 



Q[pi 



1 



2m 



(17) 



(18) 



The energy average is computed by sampling pt{x) over 
a set of trial points {xi\ and then moving the trial points 
along the conjugate gradient of 



E{x)^Vi{x) + Q[pt]{x). 



(19) 



After each conjugate gradient step, a new set of Cn coeffi- 
cients are determined according to an expectation maxi- 
mization criteria such that the new trial density provides 
the best A^-Gaussian approximation to the actual prob- 
ability distribution function sampled by the new set of 
trial points. The procedure is repeated until 5{E) = 0. 
In doing so, we simultaneously minimize the energy 
and optimize the trial function. Since the ground state 
is assumed to be node-less, we will not encounter the 
singularities and numerical instabilities associated with 
other Bohmian equations of motion based approaches. 
[H Hi, m Hi, [H [13 Moreover, the approach has been 
extended to very high-dimensions and to finite temper- 
ature by Derrickson and Bittncr in their studies of the 
structure and thermodynamics of rare gas clusters with 
up to 130 atoms. [2l|, ^ 
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IV. TEST CASE: TUNNELING IN A DOUBLE 
WELL POTENTIAL 

As a non-trivial test case, consider the tunneling of a 
particle between two minima of a symmetric double po- 
tential well. One can estimate the tunneling splitting us- 
ing semi-classical techniques by assuming that the ground 
and excited states are given by the approximate form 



V (cm"') 



1p± = -^{<j)o{x) ± (f>o{-x)) 



(20) 



where (j)o is the lowest energy state in the right-hand well 
in the limit the wells are infinitely far apart. If we assume 
the localized states {(po) to be gaussian, then 



and we can write the superpotential as 



W 



h[3 [x — Xo iSL,'ah{2xXol3)) . 



(21) 



(22) 



From this, one can easily determine both the original 
potential and the partner potential as 



Vi,: 



m 



-.W 



(2{x — Xo tanh(2xa;o/3))^ 



± {2xlsech'^i2xXo(3)-l) 



(23) 



(24) 



While the Vi potential has the characteristic double min- 
ima giving rise to a tunneling doublet, the SUSY partner 
potential V2 has a central dimple which in the limit of 
Xo 00 becomes a (5-function which produces an un- 
paired and node-less ground state. [Tj] Using Eq. [TTl 



',(1) 



(X A^ipl^"^ which now has a single 



one obtains tpl 
node at a; = 0. 

For a computational example, we take the double well 
potential to be of the form 



Viix) 



bx^ +Eo. 



(25) 



with a = 438.9cm-V(6o/ir2), b ^ 877.8cm- V(6o/ir)4, 
and Eo = — ISl.lcm^^ which (for m — uih ) gives 
rise to exactly two states at below the barrier separat- 
ing the two minima with a tunneling splitting of 59.32 
cm"^ as computed using a discrete variable representa- 
tion (DVR) approach. ^23] For the calculations reported 
here, we used Up = 1000 sample points and TV = 15 
Gaussians and in the expansion of pt{x) to converge the 
ground state. This converged the ground state to 1 : 10~^ 
in terms of the energy. This, in itself is encouraging since 
the accuracy of typical Monte Carlo evaluation would be 
~ 3% in terms of the energy. Plots of Vi and the 

converged ground state is shown in Fig. [TJ 

The partner potential V2 = + hW'/V2m, can 
be constructed once we know the superpotential, W{x). 
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FIG. 1: (a) Model double well potential(blue) and partner 
potential (purple). The energies of the tunneling doublets 
are indicated by the horizontal lines a,t V — Ocm~^ and 
V = 59.32 cm~^ indicating the positions of the sub-barrier 
tunneling doublet, (b) Final ground state density (blue) su- 
perimposed over the Gaussians used in its expansion, (purple) 



Here, we require an accurate evaluation of the ground 
state density and its first two log-derivatives. The advan- 
tage of our computational scheme is that one can eval- 
uate these analytically for a given set of coefficients. In 
Fig. [T^ we show the partner potential derived from the 
ground-state density. Where as the original Vi poten- 
tial exhibits the double well structure with minima near 
Xo = ±1 , the V2 partner potential has a pronounced dip 
about x = 0. Consequently, its ground-state should have 
a simple "gaussian" -like form peaked about the origin. 

Once we determined an accurate representation of 
the partner potential, it is now a trivial matter to re- 
introduce the partner potential into the optimization 
routes. The ground state converges easily and is shown in 
Fig. [2^ along with its gaussians. After 1000 CG steps, the 
converged energy is within 0.1% of the exact tunneling 
splitting for this model system. Again, this is an order of 
magnitude better than the Xj ^Jn^ error associated with 
a simple Monte Carlo sampling. Furthermore, Fig. [2)3 



(2) 



shows i/'i oc A^Vo computed using the converged pj,' 
density. As anticipated, it shows the proper symmetry 
and nodal position. 

By symmetry, one expects the node to lie precisely 

at the origin. However, since we have not imposed any 

symmetry restriction or bias on our numerical method, 

the position of the node provides a sensitive test of the 

(2) 

convergence of the trial density for pQ . In the exam- 
ple shown in Figl^l the location of the node oscillates 
about the origin and appears to converge exponentially 
with number of CG steps. This is remarkably good con- 
sidering that this is ultimately determined by the qual- 
ity of the 3rd and 4th derivatives of ph^ since these ap- 
pear when computing the conjugate gradient of Vi- We 
have tested this approach on a number of other one- 
dimensional bound-state problems with similar success. 
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FIG. 2: (a) Ground state density of the partner Hamiltonian 
H2 (blue) superimposed over its individual Gaussian compo- 
nents, (b) Excited state derived from the ground state 



of the partner potential, ipo 

node position 



(2) 
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FIG. 3: Location of excited state node for the last 600 CG 
steps. 



V. EXTENSION TO HIGHER DIMENSIONS 

Having demonstrated that the SUSY approach can be 
used to compute excitation energies and wave functions 
starting from a Monte Carlo approach, the immediate 
next step is to extend this to arbitrarily higher dimen- 
sions. To move beyond one dimensional SUSY, loffe and 
coworkers have explored the use of higher-order charge 
operators [H, [2^ [5^, [23] , and Kravchenko has explored 
the use of Clifford algebras [28!| . Unfortunately, this is 
difficult to do in general. The reason being that the 
Riccati factorization of the one-dimensional Schrodinger 
equation does not extend easily to higher dimensions. 
One remedy is write the charge operators as vectors 
q = {+d + W) and with ^ {^d + W f as the ad- 
joint charge operator. The original Schrodinger operator 
is then constructed as an inner-product 



Working through the vector product produces the 
Schrodinger equation 

Hi(j)= {-V^ + - {V ■]¥))(!) ^0 (27) 

and a Riccati equation of the form 

U{x) = VK^ - divW. (28) 

For a 2d harmonic oscillator, we would obtain a vector 
superpotential of the form 

W = Jiygrad^'^ = {x, y) = {W,, Wy) (29) 

Let us look at the divM^ part closer: If we us the form 
that W = -VlnV', then -V • Vlni/' = -V^lni/' which 
for the 2D oscillator results in divW^ = 2. Thus, 



- AbfW = (x^ +r) - 2 



(30) 



which agrees with the original symmetric harmonic po- 
tential. Now, we write the scaled partner potential as 



U2 = W^ + dwW = (x^ + y2) + 2. 



(31) 



This is equivalent to the original potential shifted by a 
constant amount. 



U2 = Ui+ 4. 



(32) 



(26) 



The ground state in this potential would be have the 
same energy as the states of the original potential with 
quantum numbers n + m = 2. Consequently, even with 
the this nai've factorization, one can in principle obtain 
excitation energies for higher dimensional systems, but 
there is no assurance that one can reproduce the entire 
spectrum of states. 

The problem is neither Hamiltonian H2 nor its associ- 
ated potential U2 is given correctly by the form implied 
by Eq. [27] and Eq. [211 Rather, the correct approach 
is to write the H2 Hamiltonian as a tensor by taking 
the outer product of the charges H2 = qq^ rather than 
as a scaler q-q^. At first this seems unwieldy and un- 
likely to lead anywhere since the wave function solutions 
in the second sector are now vectors rather than scalers. 
However, rather than adding an undue complexity to the 
problem, it actually simplifies matters considerably. As 
we demonstrate in a forthcoming paper, this tensor fac- 
torization preserves the SUSY algebraic structure and 
produces excitation energies for any dimensional SUSY 
system. Moreover, this produces a scaleri— > tensor 1-^ 
scaler hierarchy as one moves to higher excitations. [1^ 



VI. DISCUSSION 

In brief, we have used the ideas of SUSY quan- 
tum mechanics to obtain excitation energies and excited 
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state wavefunctions within the context of a Monte Carlo 
scheme. This was accomphshed without pre-specifying 
the location of nodes or restriction to a specific symme- 
try. While it is clear that one could continue to determine 
the complete spectrum of Hi , the real challenge is to ex- 
tend this technique to higher dimensions. Furthermore, 
the extension to multi-Fermion systems may be accom- 
plished through the use of the Gaussian Monte Carlo 
method in which any quantum state can be expressed as 
a real probability distribution, [s^, HH We offer this pa- 
per as the starting point for stimulating interest in devel- 
oping numerical techniques based upon SUSY quantum 



mechanics. 
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